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Abstract 



This is the first in a series of papers describing a 3 + 1 computational scheme 

\ for the numerical simulation of dynamic spherically symmetric black hole 

spacetimes. In this paper we discuss the construction of dynamic black hole 

qq ' initial data slices using York's conformal-decomposition algorithm in its most 

general form, where no restrictions are placed on K (the trace of the extrinsic 

curvature) and hence the full 4- vector nonlinear York equations must be solved 

00 ! numerically. 
<3\ 



To construct an initial data slice, we begin with a known black hole slice 
(e.g. a slice of Schwarzschild or Kerr spacetime), perturb this via some Ansatz 
(e.g. the addition of a suitable Gaussian to one of the coordinate components 
bJ[), of the 3-metric, extrinsic curvature, or matter field variables), apply the York 

decomposition (using a further Ansatz for the inner boundary conditions) 
to project the perturbed field variables back into the constraint hypersur- 



$— i ' face, and finally optionally apply a numerical 3-coordinate transformation 

to restore any desired form for the spatial coordinates (e.g. an areal radial 
coordinate). 

In comparison to other initial data algorithms, the key advantage of this 
algorithm is its flexibility: K is unrestricted, allowing the use of whatever slic- 
ing is most suitable for (say) a time evolution. This algorithm also offers great 
flexibility in controlling the physical content of the initial data, while placing 
no restrictions on the type of matter fields, or on spacetime's symmetries or 
lack thereof. 

We have implemented this algorithm for the spherically symmetric scalar 
field system. We present numerical results for a number of asymptotically 
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flat Eddington-Finkelstein-like initial data slices containing black holes sur- 
rounded by scalar field shells, the latter with masses ranging from as low as 
0.17 to as high as 17 times the black hole mass. In all cases we find that 
the computed slices are very accurate: Using 4th order finite differencing on 
smoothly nonuniform grids with resolutions of Ar/r m 0.02 (0.01) near the 
perturbations, and Gaussian perturbations yielding scalar field shells with rel- 
ative width a/r ~ | and with about 2/3 the black hole mass, the numerically 
computed energy and momentum constraints for the final slices are ^$ 10 -8 
(10 -9 ) and < 10~ 9 (10 -10 ) in magnitude respectively. For similarly perturbed 
numerically computed "warped" vacuum slices, the Misner-Sharp mass func- 
tion is independent of position, and the 4-Riemann quadratic curvature invari- 
ant R a bcdR abcd is equal to its (position-dependent) Schwarzschild-spacetime 
value, to within relative errors of <10" 5 (5xl0~ 7 ) and <3xl0~ 4 (5xl0~ 5 ) 
respectively. All these errors show the expected 0((Ar) 4 ) scaling with grid 
resolution Ar. 

Finally, we briefly discuss the errors incurred when interpolating data from 
one grid to another, as in numerical coordinate transformation or horizon 
finding. We show that for the usual moving-local-interpolation schemes, even 
for smooth functions the interpolation error is not smooth. If interpolated data 
is to be differentiated, we argue that either the interpolation order should 
be raised to compensate for the non-smoothness, or explicitly smoothness- 
preserving interpolation schemes should be used. 

04.25.Dm, 02.70.Bf, 02.60.Lj, 02.60.Ed, 04.40.Nr 
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I. INTRODUCTION 



Dynamic black hole spacetimes are interesting both astrophysically (e.g. in the study 
of supernova core collapse, active-galactic-nuclei engines, and binary neutron-star or black- 
hole decay and coalescence) and as laboratories in which to study the physics of strong- 
field gravitation itself (e.g. horizons, causality, cosmic censorship, and self-similarity). Such 
systems are time-dependent, strongly relativistic, and often highly asymmetric, and so are 
hard to study in detail except by numerical methods. 

This is the first in a series of papers describing a 3 + 1 computational scheme for the 
numerical simulation of dynamic black hole spacetimes. We discuss both our general com- 
putational scheme (which should be applicable to a wide range of dynamic black hole space- 
times), and its application to a specific test system, that of an asymptotically flat spherically 
symmetric spacetime containing a black hole surrounded by a (massless, minimally coupled) 
scalar field. 1 Accordingly, we first work out our formalism and describe our algorithms in 
generic terms (making no assumptions about spacetime symmetries or what types of matter 
fields may be present), and only then specialize first to the scalar field system, and finally 
to the spherically symmetric scalar field system. 

In this paper we focus on the "initial data problem" of constructing suitable initial data 
on a Cauchy surface. 2 That is, we consider the problem of computing the 3-metric and 
extrinsic curvature for an initial data slice, such that the constraints are satisfied and the 
slice has a desired physical content. 

For black hole spacetimes, we have the additional problem of dealing with singularities. 
Traditionally this has been done using freezing slicings (see, for example, Smarr (1984); 
Shapiro and Teukolsky (1986)), but more recently attention has focused on the "black hole 
exclusion" or "horizon boundary condition" technique 3 (see, for example, Thornburg (1985, 
1987, 1991); Seidel and Suen (1992); Thornburg (1993); Scheel et al. (1995a); Anninos 
et al. (1995); Scheel et al. (1995b); Marsa (1995); Marsa and Choptuik (1996); Scheel et al. 
(1997)), where the slices penetrate the event horizon and meet the singularity or singularities, 
but in each slice the computational domain (numerical grid) excludes or omits roughly the 



1 We have previously (Thornburg (1993)) applied a similar computational scheme to a different 
test system, that of an asymptotically fiat axisymmetric vacuum spacetime containing a black hole 
surrounded by gravitational radiation. Unfortunately, for that system we have been unable to 
obtain time evolutions free of finite differencing instabilities. 

2 The term "initial value problem" is sometimes used here. We find this terminology undesirable 
because it conflicts with standard usage elsewhere (e.g. in the study of differential equations, and 
indeed in this journal's own Physics and Astronomy Classification Scheme), where it refers to the 
evolution of initial data, i.e. the time integration of a set of ODEs or PDEs given values of the 
state variables at some initial time. 

3 We prefer the former term as focusing attention on the underlying process - the exclusion of 
the black hole from the computational domain - rather than on the particular boundary condition 
used to implement it. 
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black hole region. (The precise region excluded varies between different black-hole-exclusion 
computational schemes.) 

The black hole exclusion technique requires special consideration when constructing ini- 
tial data, because the boundary of the excluded region becomes a new "inner" boundary 
for the computational domain, typically on or just inside the (each) black hole's outermost 
apparent horizon. 4 This new boundary requires corresponding boundary conditions for the 
(elliptic) constraint equations. 

Because the computational domain includes the (event) horizon, the black hole exclusion 
technique requires that the slicing and the 3 + 1 field variables all be nonsingular and smooth 
there, and in fact throughout some neighbourhood of the horizon. This generally rules out 
maximal slices such as the Schwarzschild (Boyer-Lindquist) slices in Schwarzschild (Kerr) 
spacetimes; Eddington-Finkelstein (Kerr) or similar slices are usually used instead. For our 
purposes, the key property of these latter slices is that they have K, the trace of the extrinsic 
curvature, nonzero and spatially variable over most or all of the slices. 

Our slices having generic K has an important impact on the initial data construction 
itself: York's conformal-decomposition technique is widely used for solving the initial data 
problem, and we use it here, but historically it's usually been used in one of several special 
cases in which the York equations simplify greatly. However, as discussed in section IV, 
none of these special cases apply for slices with generic K, so we must (numerically) solve 
the York equations in their full form. 

In the remainder of this paper, we first summarize our notation (section II), then dis- 
cuss the York conformal decomposition, its various special cases, our general methods for 
solving its equations, and its outer boundary conditions (sections III- VI). We then describe 
our Ansatze for choosing the inputs to and the inner boundary conditions for the York de- 
composition (section VII). We next describe our diagnostics for assessing the accuracy of 
our initial data and studying its physical content (section VIII), introduce our spherically- 
symmetric-scalar- field test system (section IX), describe our numerical methods (section X), 
and present some sample results from a numerical code incorporating these techniques (sec- 
tion XI). We end the main body of the paper with some conclusions and directions for 
further research (section XII). In the appendices we tabulate some of the 3 + 1 field vari- 
ables for an Eddington-Finkelstein slice of Schwarzschild spacetime (appendix A), tabulate 
the detailed equations for our spherically-symmetric-scalar- field test system (appendix B), 
discuss the derivation of our 3 + 1 form of the Misner-Sharp mass function (appendix C), de- 
scribe our numerical-coordinate-transformation algorithm (appendix D), describe our meth- 
ods for quantitatively testing finite differencing convergence (appendix E), and discuss the 
non-smoothness of errors when interpolating data from one grid to another (appendix F). 



4 As discussed elsewhere (Thornburg (1993, 1998)), our particular form of the black hole exclusion 
technique currently places the inner boundary considerably inside the apparent horizon (typically 
at ~ 75% of the horizon radius), but the details of this are unimportant here. 
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II. NOTATION 



We generally follow the sign and notation conventions of Misner et al. (1973), with a 
(— , +, +, +) spacetime metric signature, G — c — 1 units, and all masses and coordinate 
distances also taken as dimensionless. We assume the usual Einstein summation convention 
for all repeated indices, and we use the Penrose abstract-index notation, as described by 
(e.g.) Wald (1984). However, for pedagogical convenience we often blur the distinction 
between a tensor as an abstract geometrical object and the vector or matrix of a tensor's 
coordinate components. 

We use the standard 3 + 1 formalism of Arnowitt et al. (1962) (see York (1979, 1983) for 
recent reviews). For our spherically-symmetric-scalar-field test system, we use coordinates 
(t, r, 9, 0), with 9 and <j) the usual spherical-symmetry coordinates. However, we leave t and 
r arbitrary, i.e. we make no assumptions about the choice of the lapse function or (the radial 
component of) the shift vector. 

The distinction between 3- and 4-tensors is usually clear from context, but where am- 
biguity might arise we use prefixes ® and ( 4 ) respectively, as in and ^R. Any tensor 
without a prefix is by default a 3-tensor. £ v denotes the Lie derivative operator along the 
4- vector field v a . 5 % j is the usual Kronecker delta. 

We use abed for spacetime (4-) indices, and d a denotes the spacetime coordinate par- 
tial derivative operator d/dx a . g a t denotes the spacetime metric and V a the associated 
4-covariant derivative operator. 

We use ijkl for spatial (3-) indices. di denotes the spatial coordinate partial derivative 
operator d/dx l . gij denotes the 3-metric of a slice, V« the associated 3-covariant derivative 
operator, and g the determinant of the matrix of g^s coordinate components, a and f3 l 
denote the 3 + 1 lapse function and shift vector respectively. n a denotes the (timelike) 
future pointing unit normal to the slices. = \£ n gij = —ViTij denotes the 3-extrinsic 
curvature of a slice, and K = its trace, p = n a n b T a b and j % = —n a T at denote the 
locally measured energy and 3-momentum densities respectively. denotes the spatial 
stress-energy tensor, and T = its trace. 

The 3 + 1 energy and momentum constraints are thus 

C=(R- KijK ij + K 2 ) - (l6rrp) = (la) 
C l = (VjK ij - V l K) - (87Tj*) =0 (lb) 

respectively, where R denotes the 3-Ricci scalar computed from g^ in the usual manner. 

diag[- ■ •] denotes the diagonal matrix with the specified diagonal elements. 
Gaussian(x=A, a—B) denotes the Gaussian function exp(— \z 2 ) ) where z = (x — A)/B. 

III. THE YORK DECOMPOSITION 

York (1971, 1972, 1973); 6 Murchadha and York (1974a, 1974b) have given an elegant 
theoretical analysis of the 3 + 1 initial data problem, which both clarifies its physical meaning 
and supplies a practical algorithm for its solution. (Isenberg (1979) also discusses several re- 
lated topics related to the canonical formulation of gravity.) This "conformal decomposition" 
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approach to the initial data problem is reviewed by York (1979, 1983, 1985); York and Piran 
(1982). Choptuik (1982) also gives a very clear presentation of the York decomposition. 

For present purposes, we view the York decomposition as defining a projection operator 
y which projects the field variables into the constraint hypersurface. In detail, y is defined 
as follows: 

Given an initial ("base") 3-metric gij, extrinsic curvature Kij, and energy and 3- 
momentum densities p and j l , we first split into its determinant g and the confor- 
mal metric hij = g~^ 3 gij, and into its trace K and the tracefree extrinsic curvature 
Eij = — \gijK. We then define the covariant derivative Vj and compute the 3-Ricci 
tensor and scalar i?^ and R, all using the base metric g^ in the usual manner. 

Next, we solve the nonlinear elliptic "York equation" 

ViV** = + ±K 2 ¥ J - iFijF^y- 7 - 2npty- 3 (2a) 
(A e QY = l{V l K)^ 6 - VjE lj + 87rf (2b) 

for the 4- vector 5 Y = (*,fi < ), where F ij = E ij + (£tt) ij , and where the York "vector 
divergence" and "vector Laplacian" operators £ and are defined by 

{tXf = V l X j + V j X l - \g ij V k X k (3a) 

(A e xy = Wjiexf (3b) 

= Vj V j X i + lV i V j X j + R!jX j (3c) 

for any (arbitrary) contravariant 3- vector X 1 . We discuss boundary conditions for the York 
equation in sections VI and VII. 

Finally, we define the output field variables (g^, K^, p, f) out = y(gij,K~ij,p,j z ) by 



\9ij) out 




(4a) 


(^)out 


= ^~ 2 Fij + \K^ 9lJ 


(4b) 


(P)out 


= ^- 8 p 


(4c) 


U )out 




(4d) 



Although it's not obvious from the form of the equations, the output field variables 
determined in this way do indeed satisfy the constraints. 

IV. SPECIAL CASES OF THE YORK DECOMPOSITION 

As discussed by York (1979); Bowen and York (1980); York and Piran (1982); York (1983, 
1989) and reviewed by Thornburg (1993), section A5.1, the York decomposition has several 
important special cases in which the general York equations (2) simplify considerably: 



5 Here and throughout our discussion of the York decomposition and equations, we use the term 
"4- vector" only in the sense of "4-tuple of real- valued functions", not in the sense of "rank 1 
4-tensor" . 
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• If the base field variables define a maximal slice (K = 0), or more generally if K is 
spatially constant, then the 4-vector nonlinear elliptic system (2) decouples into sepa- 
rate linear 3-vector (momentum) and nonlinear scalar (energy) constraint equations, 
which may be solved sequentially. 

• If the base field variables define a slice which is both maximal and 3-conformally flat 
(gij = § 4 fij for some conformal factor <3> and flat 3-metric fy), the York equations 
simplify still further. In particular, for this case Bowen (1979b, 1979a); Bowen and 
York (1980); Bowen (1982) have found the (analytical) general solution for the (linear) 
momentum constraint equation (2b) subject to outer boundary conditions appropriate 
for an asymptotically flat spacetime. This means that only the (nonlinear) energy 
constraint equation (2a) need be solved numerically. Bowen (op. cit.) has also found 
several Ansatze for choosing the base extrinsic curvature in this case, such that the 
final field variables represent an initial slice containing a single black hole with freely 
specifiable mass, momentum, and spin. Since the momentum constraint equation (2b) 
is linear here, these base may be superimposed to obtain base values for 
multiple-black-hole slices. We discuss this case further below. 

• If the base field variables define a slice which is time-symmetric (K^ = 0) and 3- 
conformally flat, then the full 4-vector York equation (2) has a fully-analytical asymp- 
totically flat solution, first explicitly found by Misner (1960), following an earlier sug- 
gestion of Einstein and Rosen (1935). The Misner solution, first studied numerically by 
Brill and Lindquist (1963), represents an arbitrary number of momentarily stationary 
"generalized Schwarzschild" black holes, each with freely specifiable "mass". 7 

The special case where the initial slice is maximal and 3-conformally flat has been used 
in most past numerical work on the initial data problem. In particular, several authors, 
notably York and Piran (1982); Choptuik (1982); Dubai (1992), have used this technique to 
construct initial data slices containing single black holes with varying spins, surrounded by 
various amounts of gravitational radiation. 

There have also been several extensions of this work to the multiple black hole case. 
Kulkarni et al. (1983); Kulkarni (1984b); York (1984); Bowen (1985) have used infinite 
series of "image charges" to construct a semianalytical solution for maximal 3-conformally 
flat base field variables satisfying "reflection symmetric" or "inversion symmetric" inner 
boundary conditions on the surfaces of each of N black holes, and Bowen and York (1980); 
Bowen et al. (1984); Kulkarni (1984a) have specialized this approach to the 2-black-hole case. 
Using these latter solutions, various authors, for example Rauber (1985, 1986); Cook (1990, 
1991); Dubai et al. (1992); Cook et al. (1993); Matzner et al. (1998), have constructed a 
variety of numerical solutions for \I/ corresponding to initial data slices containing two black 
holes, each with freely specifiable 6 "mass", "momentum", and "spin". 7 Thornburg (1985, 



6 The question of just what can and can't be freely specified for these initial slices is somewhat 
subtle. York (1989) gives a brief discussion of this; see also Matzner et al. (1998). 

7 "Mass" , "momentum" , and "spin" are in quotes here because (outside of spherical symmetry) 
these quantities aren't uniquely defined for the individual black holes, unless they're far enough 
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1987) has also constructed maximal 3-conformally-flat 2-black-hole numerical solutions for 
but using an apparent horizon boundary condition (a version of the black hole exclusion 

technique) instead of reflection symmetry. 

Oohara and Nakamura (1989) have used a modified version of Bowen's analytical solution 

for Q l to construct initial data slices modeling a binary neutron star system. Because no 

black holes are present in this case, they avoided the problem of inner boundary conditions 

altogether. 

As noted in the introduction, however, in this paper our concern is with initial slices 
where K is nonzero and spatially variable, so none of the special cases discussed in this 
section apply, i.e. we must solve the full 4-vector York equation (2) numerically. 



V. SOLVING THE YORK EQUATIONS 

The York equation (2) is a nonlinear elliptic PDE. To solve it, we first use a global 
Newton-Kantorovich "outer" iteration (Boyd (1989), appendix C), then use standard finite 
differencing methods to numerically solve each of the resulting sequence of "inner" linear 
elliptic PDEs. 

In detail, we first rewrite the York equation (2) in the form 

G° = ViV*tf - lm - ±K 2 m 5 + \F i:j F ij ^- 7 + 2tt P ^- 3 = (5a) 
G l = {A £ Q) 1 - f (V l A^ 6 + VjE ij - 8nf = (5b) 

where we define G = Q (Y) to be the left hand side 4-vector. We then linearize the continuum 
differential operator Q about the current continuum solution estimate Y, 

g(Y + SY) = Q(Y) + j[g(Y)] (5Y) + 0(||<5Y|| 2 ) (6) 

where 5Y is a finite perturbation in Y and the linear differential operator J[(?(Y)] is the 
linearization of the differential operator Q about the point Y. We then neglect the higher- 
order (nonlinear) terms in (6) and solve for the perturbation SY such that Q(Y + 5Y) = 0. 
This gives the linear elliptic PDE 

J[0(Y)]O$Y) = -G (7) 

to be solved for 8Y. Finally, we update the approximate solution via Y — > Y + 5Y, and 
repeat the iteration until ||G|| is small. 

It's convenient to rewrite the definition (6) of the Jacobian operator J in the differential 
form dG = J(dY) to more clearly express J's physical meaning of mapping infinitesimal 
changes in Y to infinitesimal changes in G. With this interpretation in mind, it's easy to see 
from (5) that the Jacobian is given by 

dG° = ViV* -\Rd&- ^K 2 ^ d^ 

- iFijF'V-* d^ + iFi^iidtif - 6vrp^- 4 d^> (8a) 
dG i = -4(V^)^ 5 dm + (A e dtty (8b) 



apart to have separate approximate asymptotically flat regions. 
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Given a reasonable numerical solution of the updating equation (7) (discussed in sec- 
tion X), we find in practice that the Newton-Kantorovich iteration converges rapidly and 
robustly, typically reducing 1 1 G 1 1 ^ to negligible levels (e.g. < 1CT 10 ) in 3-5 iterations. 

VI. OUTER BOUNDARY CONDITIONS 

The York equation (2) is an elliptic PDE, and as such needs boundary conditions. As- 
suming the black hole exclusion technique, there are N + 1 boundaries for a slice with N 
black holes: an inner boundary near each black hole's horizon, and an outer boundary near 
spatial infinity. We discuss the inner boundary conditions in section VII. 

The outer boundary conditions are determined by the desired asymptotic forms of \l/ and 
Q l at spatial infinity. For an asymptotically flat spacetime, we have 

^ = 1 (at spatial infinity) (9a) 
Q l = (at spatial infinity) (9b) 

Purely for the initial data computation (with no time evolution involved), one can treat 
the unbounded problem domain by compactifying, e.g. by introducing a new radial coordi- 
nate s = 1 — r min /r. This has been done by, for example, Choptuik (1982); Rauber (1985, 
1986). 

However, if the initial data computation is being done in conjunction with a 3 + 1 time 
evolution code (as is our main interest here), it's typically more convenient to not compactify 
the domain, but instead apply an approximate outer boundary condition at some finite outer 
boundary radius r max . When doing this, we can exploit the known asymptotic falloff of ^ 
and Q l to obtain more accurate results than would be obtained by simply applying the 
spatial-infinity condition (9) at r = r max . As discussed by York (1979, 1980); York and 
Piran (1982); York (1983); Oohara and Nakamura (1989); York (1989); Cook (1990, 1991); 
O Murchadha (1992), this gives the Robin outer boundary conditions 

(V fc ^)n fc + — = O « (at r = r max ) (10a) 

nj (5\ - \n l n k ) (tttf 1 + A (5\ - \rfn k ) Q k =oQ»0 (at r = r max ) (10b) 

[In practice, one can often still obtain adequate accuracy by replacing the vector Robin 
condition (10b) by the simple (and tensorially incorrect) application of a scalar Robin con- 
dition (analogous to (10a)) separately to each coordinate component of Q\ 

(V fc fi> fc + y = 0(^) ~0 (at r = r max ) (11) 

We use this approximation in our numerical code, and find that it gives good results.] 

[It's also useful to note that for many coordinate conditions, the asymptotic forms of ^ 
and Q l are identical in their leading terms to those for the 3 + 1 lapse function a and shift 
vector f3 l (respectively). Hence for such coordinate conditions the boundary conditions for 
\l> and Vt l are also identical to those for a and f3 l . In the context of a 3 + 1 code incorporating 
both initial data computation and time evolution, this may allow the same outer boundary 
condition code to be used for the York equation as for the coordinate conditions.] 
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VII. OUR INITIAL DATA ALGORITHM AND ANSATZE 



Given the York projection operator, our overall initial data algorithm is as follows: 

1. We begin with an exact (analytically known) black hole initial data slice, chosen ac- 
cording to some Ansatz. As discussed in sections I and IV, this will typically have K 
nonzero and spatially variable throughout the slice. 

2. We apply an arbitrary (in general constraint- violating) perturbation, again chosen 
according to some Ansatz. 

3. We apply the York decomposition to project the perturbed field variables back into 
the constraint hypersurface, solving the full 4-vector York equations numerically. Here 
the outer boundary conditions are dictated by the physics (cf. section VI), but an 
(another) Ansatz is needed to supply the inner boundary conditions. 

4. Optionally, we apply a numerical 3-coordinate transformation to restore the spatial 
coordinates within the slice to some desired form, e.g. to make the radial coordinate 
areal (g$o = r 2 ). We discuss numerical coordinate transformations in appendix D. 

In comparison to other initial data algorithms, the key advantage of using the full York 
decomposition is its flexibility: Since no restriction is placed on K, almost any slicing may 
be used. That is, with this algorithm the slicing may be chosen based on its suitability for 
(say) a black-hole-excluded time evolution, rather than being limited by the initial data 
solver. As well, the physical content of the initial slice may be controlled over a wide range 
by varying the choices for the different Ansatze. This algorithm also places no restrictions 
on what type(s) of matter field(s) may be present, and (given suitable programming for the 
numerical solution of the York equations) it works equally well in spherical symmetry, in 
axisymmetry, or in fully generic spacetimes with no continuous symmetries. 

To actually use this algorithm to construct initial data slices, we must make suitable 
choices for the various Ansatze. We now consider these choices: 

As noted in the introduction, the black hole exclusion technique requires a slicing 
which penetrates the horizon, with the slicing and the 3 + 1 field variables all nonsingular 
and smooth throughout some neighbourhood of the horizon. We thus use an Eddington- 
Finkelstein slice of Schwarzschild spacetime, or a Kerr slice of Kerr spacetime, 8 for step 1 of 
our algorithm. 

The choice of the perturbation Ansatz for step 2 of our algorithm is much less clear- 
cut, and we have experimented with a number of different choices for what field variable to 
perturb and what perturbation to apply. Provided the perturbation isn't too large, we have 



8 Recall that these slices are both defined (Misner et al. (1973), boxes 31.2 and 33.2) by taking 
the usual areal radial coordinate r, and defining t such that t + r is an ingoing null coordinate. 
Both slices are nonsingular near to and on the horizon. For reference, in appendix A we tabulate 
some of the 3+1 field variables for an Eddington-Finkelstein slice of Schwarzschild spacetime. 



10 



found that in practice essentially any smooth localized perturbation 9 in the field variables 
"works", in the sense that the York equations remain numerically well-behaved and yield 
a properly-constrained initial data set (still) containing a black hole. However, the precise 
relationship between the form of the perturbation and the nature of the resulting initial data 
slice is complicated and not (yet) easily predicted a priori. In section XI we present numerical 
results for several examples of this relationship (different choices of the perturbation), but 
much more research is needed to explore and map this relationship's phenomenology. 

To motivate the choice of an inner boundary condition Ansatz for the York equation (2) 
in step 3 of our algorithm, we first observe that if the base field variables already satisfy the 
constraints, then the York projection operator is the identity operator, i.e. the solution of 
the York equation (2) is ^ = 1 and Q l = 0. We then consider the slightly more general case 
where the base field variables are a small perturbation off the constraint hypersurface. To 
the extent that the perturbation is indeed small, or at least that its effects are small near 
the inner boundary, then \l/ = 1 and Q 1 = should still be an approximate solution of the 
York equation at the inner boundary. This suggests taking 

^ = 1 (at the inner boundary) (12a) 
Q l = (at the inner boundary) (12b) 

as our inner boundary condition Ansatz for the York equation (2). 

It's not a priori obvious that (i) this boundary condition is consistent with any nontrivial 
solutions of the York equation (2), nor that (ii) the resulting slices will (still) contain black 
holes. However, from our numerical results in section XI, both (i) and (ii) do in fact hold in 
practice. 

Moreover, although our motivation for this boundary condition took the perturbation to 
either be small, or at least to be small in its effects near the inner boundary, our numerical 
results in section XI D also show that this restriction isn't needed, i.e. that accurate results 
(initial data slices which satisfy the constraints) are still obtained even for perturbations 
which are large and/or have large effects near the inner boundary. In fact, it appears that 
essentially any reasonable inner boundary condition still yields accurate results. 

We have previously (Thornburg (1993), sections 5.3 and 9.8) suggested that the fact that 
the physics of the initial data problem doesn't specify any particular inner boundary condi- 
tion for the York equation, constitutes a serious drawback of the overall black-hole-exclusion 
technique, and that the Ansatz (12) therefore "lacks a clear physical motivation and is at 
best an ad-hoc solution" . However, we now regard this criticism as unwarranted: Physically, 
the inner boundary condition for the York equation specifies what type of black hole (or 
more generally, spacetime region, cf. our discussion of the pqw5c slices in section XI E) is 
contained within the inner boundary. This must be specified as an input to the initial data 
algorithm, and the Ansatz (12) seems quite reasonable for this purpose. 

We have made a few limited experiments with inner boundary conditions other than 
the simple Dirichlet condition (12), but much further research is needed to understand the 
relative merits of different inner boundary conditions, and their effects on the resulting initial 
data slices. 



9 We haven't yet investigated non-localized perturbations. 
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VIII. DIAGNOSTICS 



The purpose of any initial data algorithm is to construct slices which satisfy the con- 
straints, so we use the numerically computed values of the constraints (1), (i.e. their devia- 
tions C and C % from satisfaction), as diagnostics of our initial slices' accuracy. Because of 
its dependence through R on the 2nd derivatives of the 3-metric components, the energy 
constraint C is particularly sensitive (useful) for this purpose; in practice, nonzero computed 
values of C primarily measure errors in the computation of R. 

We use a number of diagnostics to study the physical content of our slices and the 
spacetimes in which they're embedded, including the coordinate position(s) of any apparent 
horizon(s) in the slice, the mass distribution in each slice (measured in spherical symmetry 
by the Misner- Sharp mass function mus or the integrated scalar field mass function m M 
defined below), the scalar field's 3-energy density p (or the radial scalar field density Ang ee p 
in spherical symmetry), K, R, the 3-Ricci quadratic curvature invariant R^R l \ the 4-Ricci 
scalar ^R, and the 4-Riemann quadratic curvature invariant R a bcdR abcd - 

As discussed by Guven and O Murchadha (1995), in spherical symmetry the Misner- 
Sharp mass function mus can be written as the volume integral of a local function p = 
p(p,j l ; gij, Kij) of the scalar field variables. We can thus define an alternative mass function 
(which should be equal to the Misner-Sharp one) by volume-integrating p, 

An g ee p^/g^dr (13) 

- min 

where mMs( r min) is the Misner-Sharp mass function at the inner grid boundary r = r m j n . 
?71ms an d m M are computed independently and in very different ways ((B18) vs. (B20) 
and (B21)), so their numerical equality at all radii is a useful consistency check and accuracy 
diagnostic. As discussed in appendix B 3, for our computational scheme m M is numerically 
better-behaved than ttims, so we generally use m M as a general-purpose mass function, and 
hereinafter we take m (with no subscripts) to denote m^. 

By contracting the spacetime Einstein equations, it's easy to see that ^R = 8n(p — T), 
so ^R plays the role of a 4-invariant diagnostic for the scalar field density. 

R a bcdR abcd is a 4-invariant diagnostic for spacetime curvature. As discussed by York 
(1989), R a bcdR abcd is given in terms of the 3 + 1 variables by 

RabcdR ahcd = 8R tj R tj - IQRiji^K 1 k K kj - KK tj ) + 8K i:j K jk K kl K H 
- lQKK ij K i k K kj + 8K 2 KijK ij - 8V {i K j]k V l K jk 

+ 167r[4Iij(iT k K kj - KK ij ) - 1 /,',// '•' 

+ (16tt) 2 [T tJ T> - \T 2 - y + \Tp] (14) 

Most of our diagnostics vary rapidly with spatial position in any black hole spacetime, 
even Schwarzschild spacetime (cf. appendix A). For example, in an Eddington-Finkelstein 
slice of a Schwarzschild spacetime, R typically varies as ~ 1/ r 4 , where r is the areal radius. To 
focus on the deviations away from a Schwarzschild / Eddington-Finkelstein slice, we usually 
normalize out the "background" Schwarzschild variation by using "relative" diagnostics, 
e.g. R/Rschw Here for any diagnostic Z, at each areal radius, -Zschw denotes the value of 
Z at that same areal radius in an Eddington-Finkelstein slice of a unit-mass Schwarzschild 
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spacetime, and more generally, for any fixed to*, Z Scllw ^ m ^ denotes the value of Z at that 
same areal radius in an Eddington-Finkelstein slice of a mass-m* Schwarzschild spacetime. 

When considering a (spherically symmetric) slice or region of a slice whose mass function 
varies significantly with position, we also use "mass-relative" diagnostics, e.g. Rj Rs c hw(m(r)) ■ 
Here for any slice S and diagnostic Z, at each areal radius r, Zschw(m(r)) denotes the value of 
Z at that same areal radius r in an Eddington-Finkelstein slice of a Schwarzschild spacetime 
of mass m(r), where m(r) is the mass function of the slice S at the areal radius r. Marsa 
(1995); Marsa and Choptuik (1996) have used the same notion of mass-relative variables in 
their initial data construction, though not with diagnostics. 



IX. SCALAR FIELD AND SPHERICAL SYMMETRY 

We have previously (Thornburg (1993), appendix 5) described the use of the same York- 
decomposition initial data problem algorithm and Ansatze described here, in constructing 
dynamic vacuum axisymmetric initial data slices containing a black hole surrounded by 
gravitational radiation. Here we discuss the algorithm and Ansatze's use in a different sys- 
tem, the construction of dynamic spherically symmetric initial data slices (each) containing 
a black hole surrounded by one or more scalar field shells. 

The spherically symmetric scalar field and similar systems have been studied by a number 
of past researchers, including (among others) analytical studies by Christoudolou (1986a, 
1986b, 1987a, 1987b, 1991, 1993), 3 + 1 studies by Choptuik (1986, 1991); Seidel and Suen 
(1992); Bernstein and Bartnik (1995); Scheel et al. (1995a, 1995b); Anninos et al. (1995); 
Marsa (1995); Marsa and Choptuik (1996), 2 + 2 studies by Goldwirth and Piran (1987); 
Goldwirth et al. (1989); Hamade and Stewart (1996), a hybrid 3 + 1 and 2 + 2 study by 
Gomez et al. (1996), and an interesting comparison of 3 + 1 and 2 + 2 methods by Choptuik 
et al. (1992). 

[However, these authors haven't used the York decomposition for constructing their ini- 
tial data: Focusing on the 3+1 studies, Choptuik (1986, 1991); Bernstein and Bartnik (1995); 
Scheel et al. (1995a, 1995b) all constructed initial data using slicing conditions which simpli- 
fied the constraints sufficiently to allow a "direct" (algebraic) solution of the constraints by 
suitably ordering the 3 + 1 equations. Using similar Eddington-Finkelstein-like slices to ours, 
Marsa (1995); Marsa and Choptuik (1996) used an Ansatz of setting K r r = (K r r) Sc hw(m(r))i 
then used a functional iteration scheme to satisfy the remaining constraint equations. Seidel 
and Suen (1992); Anninos et al. (1995) used Schwarzschild spacetime as their only gen- 
eral relativistic system, and hence needed only to use an appropriate slice of Schwarzschild 
spacetime for their initial data.] 

Turning now to our scalar field formalism, we begin with the generic case, making no 
assumptions about about any spacetime symmetries. Following Choptuik (1986, 1991); 
Marsa (1995); Marsa and Choptuik (1996), we take the scalar field to satisfy the 4-scalar 
wave equation V a V a = 0, and to have the stress-energy tensor 

4nT ab = (<9 a 0)(«9 b 0) - |^ 6 (9 C 0)(9 C 0) (15) 

We then define the 3 + 1 scalar field variables 
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Q = {d t <j>-(?V i <l>)/cL 



so that 



d t <, 
die 



aQ + P k P k 
P 



(16a) 
(16b) 



(17a) 
(17b) 



(These definitions (16) and (17) are similar, but not identical, to those of Choptuik (1986, 
1991); Marsa (1995); Marsa and Choptuik (1996).) A straightforward calculation then gives 
the 3 + 1 energy and momentum densities as 



Avp = \{P k P k + Q 2 ) 
47cji = -PiQ 

and the spatial stress-energy tensor and its trace as 

-P k P k + Q 2 ) 



\9ij( 



AnT = -\P k P k + §Q 2 
Corresponding to (4), the York-transformation scalings for p and Q are 

(^)out = 
(Q)out = *~*Q 



(18a) 
(18b) 

(19) 
(20) 



(21a) 
(21b) 



We now assume that spacetime is spherically symmetric, with the spatial coordinates 
x % = (r, 9, (j>) having the usual polar spherical topology. However, we leave r arbitrary, i.e. 
we make no assumption about the choice of (the radial component of) the shift vector. We 
take the 3 + 1 field tensors to have the coordinate components 



gij = diag 



Kij = diag 



ABB sin 2 9 
X Y Y sin 2 9 



o o 

Pi=[P 
and for the York decomposition we take 



n o o 



(22) 
(23) 
(24) 
(25) 

(26) 



Notice that we do not factor out either a conformal factor or any r 2 factors from the 
3- metric components. 10 



°Although this somewhat simplifies the 3 + 1 equations, it may degrade the accuracy of our 
computational scheme at large r. In particular, in hindsight not factoring out r 2 factors may have 
been an unwise design choice. It might be interesting to try a side-by-side accuracy comparison 
between otherwise-identical factored and unfactored schemes, but we haven't investigated this. 
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With these assumptions, it's then straightforward, though tedious, to express all the 
other 3 + 1 variables in terms of the "state variables" A, B, X, Y, P, and Q, and their 
spatial (1st and 2nd) derivatives. For reference, in appendix B we tabulate the resulting 
equations for all the 3 + 1 variables involved in our initial data algorithm. 

X. NUMERICAL METHODS 

We use standard finite differencing techniques to numerically solve the York updating 
equation (7). We describe our numerical methods in detail in a following paper (Thornburg 
(1998)); we merely summarize them here. 

For the initial data equations, we use the usual centered 5 point 4th order finite dif- 
ference molecules in the grid interior, and off-centered 5 (6) point 4th order molecules for 
1st (2nd) derivatives near the grid boundaries. We use the same finite difference molecules in 
all contexts in the equations, including "left hand side" , "right hand side" , interior equations, 
and boundary conditions. 

After finite differencing the elliptic PDE (7), we row scale the resulting linear system to 
improve its numerical conditioning, then LU decompose and solve it via LINPACK band 
matrix routines (Dongarra et al. (1979)). We use IEEE double (64 bit) precision for all 
floating point computations. 

We use a smoothly nonuniform grid, chosen to be uniform in the "warped" radial coor- 
dinate 

where r , a, b, and c are suitably chosen parameters. We don't use any sort of staggered 
grid. Note that we use w only for finite differencing - all tensor components are still taken 
with respect to the (r, 9, <fi) coordinate basis. 

By convention, we always place the inner grid boundary at w — 0, i.e. r = r = r min . All 
our results in this paper use the parameters r = 1.5, a = oo (effectively omitting the first 
term in (27)), 6 = 5, and c = 100. As shown in figure 1, for these parameters w qualitatively 
resembles a logarithmic radial coordinate in the inner part of the grid, and a uniform radial 
coordinate in the outer part of the grid. 

XI. SAMPLE RESULTS 

As shown in table I, we have computed a number of test initial data slices, using various 
combinations of numerical parameters (grid resolutions and outer boundary positions) and 
initial perturbations. Most of our discussion of these slices is independent of their numerical 
parameters, so we usually refer generically to the families of slices sharing common initial 
perturbations, and plot only the highest-resolution results within each family. 

In describing our results, we use the subscript " init " to refer to the initial slice (step 1 
in our overall initial data algorithm of section VII), " pcr t U rb" to refer to the perturbed slice 
(the result of step 2), "York" to refer to the result of the York projection operator (step 3), 
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and "final" to refer to the final result after the numerical coordinate transformation back to 
an areal radial coordinate (step 4). 

All our results presented here use an Eddington-Finkelstein slice of a unit-mass 
Schwarzschild spacetime for the initial slice (step 1), and include a final numerical 3- 
coordinate transformation to make the radial coordinate areal (step 4). 

A. Scalar Field Phenomenology 

As a first example of our initial data solver, consider the pqw5 slice, defined by the 
perturbation (step 2 in our overall initial data algorithm of section VII) 

P -> P + 0.02 x Gaussian(r init =20, cx=5) (28) 

(Q remains identically zero.) 

Figure 2 shows the York-decomposition variables * and f2 for this slice. Notice that 
because the York equation (2) is elliptic, the localized perturbation (28) results in the York 
projection operator being nontrivial (^ ^ 1, Q ^ 0) throughout almost the entire slice. 
However, for the relatively small perturbation (28), the York projection and the numeri- 
cal coordinate transformation back to an areal radial coordinate (the latter given here by 
^finai/^York = ^ 2 ) both differ from the identity by only a few percent. 

Figure 3 shows the relative 3-metric and extrinsic curvature components for this slice, and 
figure 4 shows the resulting relative 3-invariants K, R, and RijR^ . Due to the nonlocality 
of the York projection operator, all the geometric field variables deviate significantly from 
their initial (Schwarzschild / Eddington-Finkelstein) values throughout the slice. Notice also 
that there's an apparent horizon just inside r = 2, i.e. the slice does (still) contain a black 
hole. 

Figure 5 shows this slice's scalar field and mass distributions. As can be seen, the slice 
contains a shell of scalar field surrounding the black hole. The scalar field shell has a roughly 
Gaussian profile, which for the radial density 4nBp is centered at r ^ 21.8, with thickness 
a Ki 3.5. Despite the relatively small perturbation (28), the scalar field shell is fairly massive, 
approximately 0.66 times the black hole mass (approximately 40% of the slice's total mass). 
[Although it's not apparent here, time-evolving this initial data (Thornburg (1998)) shows 
that the shell is actually the superposition of two momentarily coincident shells, one ingoing 
and one outgoing. However, contrary to what one might expect, the two shells have quite 
different masses.] 

B. Spacetime and Slicing Phenomenology 

From figure 5 it's clear that away from the scalar field shell the pqw5 slice is nearly 
vacuum (all our scalar field diagnostics are very small), so by Birkhoff's theorem it must be 
nearly Schwarzschild there. More precisely, inside the scalar field shell the pqw5 slice must 
be nearly (isometric to) some slice of a Schwarzschild spacetime with mass m < = m hh ss 
0.976, while outside the scalar field shell it must be nearly (isometric to) some slice of a 
Schwarzschild spacetime with mass m> = m tota i ~ 1.617. 
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In support of this, figure 6 shows the mass-relative 4-Riemann curvature invariant 
Rabcd,R abcd for the pqw5 slice. As expected, away from the scalar field shell R a t, c dR abcd is 
almost exactly equal to its value for a Schwarzschild spacetime of matching mass (m< inside 
the shell, m > outside). 

Figure 7 shows the mass-relative 3-invariants K, R, and RijR^ for the pqw5 slice. The 
mass-relative ratios all deviate substantially from unity even in the near-vacuum regions. 
From this we conclude that although the pqw5 slice is nearly (isometric to) some slice of a 
matching-mass Schwarzschild spacetime in each of its near-vacuum regions, neither of these 
slices is our canonical Eddington-Finkelstein one. 

It might be interesting to explicitly compute the embedding of these regions of the 
slice in their respective Schwarzschild spacetimes, i.e. explicitly compute (say) each region's 
Eddington-Finkelstein or Schwarzschild time coordinate as a function of spatial position (e.g. 
areal radius r). The two time coordinates would each only be defined up to an arbitrary 
additive constant, and (being referred to different Schwarzschild spacetimes) they wouldn't 
be directly comparable, but their intra-region variation might still be interesting. However, 
we haven't as yet made such a study. 

From figures 4 and 7, it appears that at large r, K / 'Aschw 1 and R/Rschw 1, while 
J = RijRv displays the different behavior J / Jschw(m(r)) 1 and hence (cf. (All)) J / Jschw 
(^totai/^init) 2 - Examining other slices computed with a larger outer boundary radius, e.g. 
the 100o30.pqw5 and 200o30.pqw5 slices (cf. figure 10) confirms these asymptotic behaviors. 
A simple analytical argument also confirms this behavior for K: The transformation (4) in 
the York projection operator preserves K at any given event, and since VP — > 1 at large r, 
the transformation back to an areal radial coordinate preserves K J s leading order behavior 
there. Thus our initial data slices must have asymptotically the same K as the initial 
Schwarzschild / Eddington-Finkelstein slice, i.e. K / K Schw ( m . nit ) — > 1 at large r. 

Since our different 3-invariant diagnostics have different asymptotic behaviors relative to 
our canonical Eddington-Finkelstein slice, it's clear that our slice isn't (even) asymptotically 
Eddington-Finkelstein. This asymptotic "warping" of the slice may prove a complication 
for 3 + 1 time evolution calculations. It might be worth investigating whether modifications 
to our initial data algorithm (e.g. a different type of perturbation, and/or different outer 
boundary conditions) might yield slices with nicer asymptotic properties, but we haven't as 
yet done this. 

C. Accuracy 

The main error sources in a 3 + 1 code such as ours are due to the finite grid resolution 
(finite differencing truncation errors), the finite outer boundary radius, and floating point 
roundoff. As discussed in appendix E, we expect the first two of these to vary in predictable 
ways with the corresponding numerical parameters (grid resolutions and outer boundary 
positions), and we can quantitatively measure these errors by comparing computations of 
the same physical system using different numerical parameters. In the context of a finite 
differencing code such as ours, i.e. one solving elliptic PDEs such as the York equation (2), 
floating point roundoff errors are easily distinguished from other error sources by their non- 
smoothness (large and quasi-random variations from one grid point to the next). 
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As an example of our computational scheme's overall accuracy level, figure 8 shows the 
magnitude of the numerically computed energy and momentum constraints at various points 
in the initial data computation for the 100.pqw5 and 200.pqw5 slices. Consider first part (a) 
of the figure, which shows the magnitude |C| of the energy constraint, and for the moment 
focus only on the 100.pqw5 curves: |Cs c hw| shows the basic "background" accuracy level of 
our finite differencing, tailing off into floating point roundoff noise for r ^ 40. In contrast, 
|C P erturb| shows the effects of the perturbation (28): |C per turb| is large (>> |Cs c hw|) wherever 
the perturbation is significant. Finally, |C n nai| shows the effects of the York projection 
operator (and the numerical coordinate transformation back to an areal radial coordinate), 
with |Cfi na i| <C | Cperturb I wherever the perturbation is significant, and more generally |Cfi na i| 
small (comparable in magnitude to |Cschw|) everywhere in the slice. In other words, we see 
that the York projection operator does indeed act to retore the field variables back into the 
constraint hypersurface. 

To assess our computational scheme's finite differencing errors, figure 8(a) also shows 
|Cschw|j | Cperturb |, and |Cfi na i| for the 200.pqw5 slice, which (cf. table I) is computed using 
twice the grid resolution (half the grid spacing) as the 100.pqw5 slice, but is otherwise 
identical. Notice that for |Cs cnw | and |C n nai|, the 100.pqw5 and 200.pqw5 curves are both 
very similar in shape, the 200.pqw5 curve being offset downwards (smaller |C|) by a factor 
of 16 (compare with the scale bar), until both tail off into roundoff noise at large r. As 
discussed in appendix E, this constitutes strong evidence that these |C| values are indeed 
dominated by 4th order finite differencing truncation errors at small and moderate r, and 
floating point roundoff errors at large r. 11 Notice also that the overall level of the |C| values, 
i.e. of the deviation of the energy constraint from satisfaction, is very low, with |C nna i| ^ 10~ 8 
(lCT 9 ) for the 100.pqw5 (200.pqw5) slice. 

Figure 8(b) shows the magnitude of the numerically computed momentum constraint 
C r 12 at the same points in the initial data computation for these slices. (Since the pertur- 
bation (28) doesn't change C r , |Cp Crturb | = |Cg chw |, so these are plotted as the same curve.) 
From the figure it's clear that the overall level of the momentum constraints is ~ 3 orders of 
magnitude smaller than that of the energy constraint. Also, |Cg nal | is small (comparable in 
magnitude to |Cg chw |) everywhere in the slice, and both show 4th order convergence in the 
same manner as |Cschw| and |C nna i|. 

As an example of a more quantitative test of the finite differencing errors, figure 9 shows 
a scatterplot of the 200.pqw5 versus the 100.pqw5 C nna i values. As discussed in appendix E, 
for 4th order convergence all the points (except for outliers from the grid boundaries) should 
fall on a line through the origin with slope ^. As can be seen, except for the grid-boundary 
outliers, all the points do indeed fall very closely on this line. We emphasize here that 
(cf. appendix E) the plotted line is not a fit to the data, but rather an a priori prediction 



11 Due to boundary finite differencing effects (cf. appendix E), C and C r are both roughly a factor 
of 10 larger at the grid boundaries than at nearby non-boundary grid points. However, even the 
boundary C and C r values still show 4th order convergence. 

12 In spherical symmetry C r is the only nonzero component of the momentum constraint vector 
C\ cf. (B16). 
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with no adjustable parameters, making this a very stringent test. 

Errors due to the finite outer boundary radius r max are more problematical: these enter 
through the outer boundary conditions (10a) and (11), which essentially fix the asymptotic 
flatness of the initial data slice. Without a clear-cut measurement of the slice's embedding 
in the far-field Schwarzschild region (cf. section XI B), asymptotic-flatness errors are difficult 
to measure. Moreover, it's not immediately obvious just how (i.e. with what power of r max ) 
these errors should scale with r max . 

However, to give an indication of the typical magnitudes of these errors, the relative 
changes 5K/K, 5R/R, 5 J/ J, and 51 /I (where J = R ij R ij and / = R abcd R abcd ) between the 
same positions in (common to) the pqw5 slices with u> max = 4 (r max ps 248) and u> max = 10 
(r max m 813), are all on the order of 10~ 5 , while the relative changes 5m/ m in the computed 
black hole and total spacetime masses are ~3xl0~ 6 . The corresponding changes between 
the slices with u> max = 10 (r max ~ 813) and w max = 30 (r max ps 2 780), are all about a 
factor of 10 smaller, and in general the outer boundary errors appear to scale approximately 
as ~l/r max 2 . These are all continuum effects: with one exception discussed below, the 
differences between different-resolution slices with the same r max , are all much smaller than 
the differences between slices with different r max . 

The one significant non-continuum effect apparent in these slices is the presence of sig- 
nificant floating-point roundoff noise in our curvature diagnostics, particularly in R, at large 
r. As an example of this, figure 10 shows the relative and mass-relative values of various em- 
bedding and curvature diagnostics for the outer parts of the 100o30.pqw5 and 200o30.pqw5 
slices. While the K, RijR % \ and and R a bcdR ahcd values show no visible noise, the R values 
show serious noise, with the relative noise amplitude 5R/ Rs c hw roughly 4 times larger in the 
200o30.pqw5 slice, and in both slices growing rapidly with r to ~5% (20%) at the outer 
boundary (r max ps 2 780) of the 100o30.pqw5 (200o30.pqw5) slice. The rapid growth of the 
relative noise amplitude with r is actually entirely due to the rapid (~ 1/r 4 ) falloff in -Rschw^ 
within each slice the actual absolute noise level in R is essentially constant for the range of 
r in the figure. 

Comparing the computations (appendix B) of the different diagnostics plotted in fig- 
ure 10, we see that K is computed from the 3-metric and extrinsic curvature components 
by algebraic operations only, while R, RijR 1 *, and R a b c dR abcd all depend on (2nd) numer- 
ical derivatives of the 3-metric components. Numerical differentiation is well known as a 
noise-amplifying process, and will inherently amplify normal low-level roundoff noise in the 
3-metric components. To test whether this process can account for the observed noise in 
R, figure 10 also shows the relative value of R &t ± e/(Aw) 2 , where i? fit is an empirical 
least-squares fit to the R values, and the ±e/(Aw) 2 term models the effects on R of 0(e) 
roundoff errors in the 3-metric components being 2nd-differentiated, with e = 6xl0~ 19 an 
"eyeball-fitted" constant parameter, and Aw the grid spacing in our w nonuniform gridding 
coordinate (cf. section X). From the figure it's clear that this model does in fact nicely 
account for the r- and Aw-dependence of i?'s roundoff noise. 

Although this explains the noise in R, it raises the further question of why R^R^ and 
RabcdR abcd don't show similar effects. (A careful analysis does show some roundoff noise in 
RijR'i and R a b c dR abcd , but only at levels many orders of magnitude below that in R.) We're 
not certain, but we suspect the explanation lies in the detailed form of the respective equa- 
tions, with the offending 2nd derivative terms being (presumably) relatively less important 
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in the computations of RijR^ and R a bcdR abcd than in R. In any case, though, the absolute 
magnitude of even i?'s roundoff noise is still very small, ±e/(Aw) 2 ~ 6xl0" 15 (2xl0~ 14 ) 
for the 100o30.pqw5 (200o30.pqw5) slice. Although R appears as a coefficient in the York 
equation (2), the offending term is linear (so -R's noise should mainly average out), and 
the noise doesn't appear to be a serious error source in the solutions \I/ and Q\ or in the 
computed 3-metric and extrinsic curvature components. 

D. Inner Boundary Condition 

Thus far we have considered only the perturbation (28), which is relatively small at the 
inner boundary (5P inner /5P max ~ 10~ 3 ). As discussed in section VII, our motivation for the 
inner boundary condition (f2) is actually predicated on the perturbation having negligible 
effects there. 

It's thus of interest to determine whether and how well our initial data algorithm works 
for perturbations which don't satisfy this restriction, i.e. for perturbations which have non- 
negligible effects at the inner boundary. As an example of such a perturbation, consider the 
pqw5i slice. As shown in table I, this uses the perturbation 

P -> P + 0.02 x Gaussian(r init =f 0, (7=5) (29) 

instead of (28). The perturbation (29) has the same amplitude and width as (28), but 
is centered at r init = 10 instead of r init = 20, so it's fairly large at the inner boundary 

(5-Pinner/ ^-^max ~ 0.24). 

Figure II shows |Cs c hw|, |Cpcrturb|, and |Cg na i| for the pqw5i slice, analogously to fig- 
ure 8(a) for the pqw5 slice. It's clear that despite the pqw5i slice's substantial pertur- 
bation at the inner boundary, its initial data computation remains fully accurate, with 
|Cfinai| <C |C pcrtur b| wherever the perturbation is significant, and |Cfi na i| small (comparable in 
magnitude to |Cs c hw|) everywhere in the slice, including near to and at the inner boundary. 

Moreover, comparing the |Cschw| and |Cfi na i| curves in figure II, we see that the 100.pqw5i 
and 200.pqw5i curves are again both very similar in shape, the 200.pqw5i curve being offset 
downwards (smaller |C|) by a factor of 16, until both tail off into roundoff noise at large 
r. Just as for the pqw5 slice (cf. our discussion of figure 8 in section XI C), this implies 
that these errors (the nonzero values of C) are dominated by (the expected) 4th order finite 
differencing truncation errors at small and moderate r, and floating point roundoff errors 
at large r. 13 The same is also true of C r for these slices (details of the analysis omitted for 
brevity) . 

These pqw5i results demonstrate that despite our inner boundary condition (12) being 
motivated by - and only for - the special case where the perturbation's effects are negligible 
at the inner boundary, we actually obtain fully accurate initial data even for perturbations 
which are large there. 



13 The pqw5i \C\ values in figure 11 also show similar boundary finite differencing anomalies 
(cf. appendix E), to the pqw5 |C| values in figure 8 (cf. footnote 11). 
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E. Other Perturbations 

As shown in table I, we have also experimented with several alternative variants of the 
perturbation (28) for step 2 in our algorithm summary of section VII. 

Figure 12 shows the effect of increasing the perturbation amplitude from 0.02 (the pqw5 
slice discussed above) to 0.05 (the pqw5b slice) or 0.1 (the pqw5c slice). With one excep- 
tion disussed below, all three slices are qualitatively fairly similar. In going from the pqw5 
slice to the pqw5b and pqw5c slices, the main detailed differences are that as a result of 
the larger perturbations, the scalar field shells become much more massive (from approx- 
imately 0.66 to 3.9 and 17 times the black hole mass, or equivalently from approximately 
40% to 80% and 94% of the total spacetime mass, respectively), and are also somewhat 
larger and thicker. As well, the deviations (not shown here) of all the geometric field vari- 
ables from their Schwarzschild / Eddington-Finkelstein values become much larger. Despite 
these stronger nonlinear effects, we experienced no difficulties with poor convergence of the 
Newton-Kantorovich iteration for these models. 14 

The pqw5c slice shows a new effect: although the mass function still shows a significant 
mass contained within the numerical grid's inner boundary (at r = 1.5), here there's no 
longer any apparent horizon within the numerical grid. The pqw5c slice's "black hole" 
mass shown in table I thus no longer represents the mass of an actual black hole, but rather 
merely of the region within the inner boundary. Since r « 2.6m > 2m on the inner boundary 
this region may, but need not necessarily, contain a black hole. (Even with the constraints 
imposed, the continuation of the field variables into the origin is obviously non- unique.) 

Since the inner grid boundary for this slice doesn't lie within an apparent horizon, this 
numerical data is unsuitable for a black-hole-excluded time evolution. However, viewed 
purely as initial data, it remains fully accurate, i.e. the constraints are small everywhere in 
the numerical grid. 

Figure 13 shows the effects of using a narrow (width-1) but high-amplitude (0.1) per- 
turbation (the pqwl slice). This slice qualitatively resembles the pqw5b slice: it contain a 
black hole surrounded by a moderately massive scalar field shell (approximately 3.0 times 
the black hole mass, or equivalently approximately 75% of the total spacetime mass). How- 
ever, corresponding to the narrow perturbation, here the scalar field shell is very thin, with 
a fractional width (cr/r centCT for a Gaussian fit to the radial scalar field density AnBp) of only 
0.033. 

The pw5+qw3 slice shows the effect of applying perturbations to both P and Q, instead 
of P alone as in the previous slices. This slice is qualitatively fairly similar to the pqw5b 
slice, except that since both P and Q are nonzero, j r is nonzero as well as p. 

Having j r nonzero makes this slice particularly useful for testing the equality of the 
Misner-Sharp and integrated-scalar- field mass functions mMs and m M (cf. section VIII), 
since now both terms in the definition (B20) of p contribute, instead of just the p term 
as in all the previous slices. Figure 14 shows the relative deviation of m MS from m M for 



14 In retrospect this isn't too surprising, as even for the pqw5c slice the peak values of \& — 1 and 
are still only «0.22 and 0.16 respectively, so the nonlinearities are still fairly weak in an absolute 
sense. 
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this slice, for both grid resolutions. Both deviations are very small across their slices, with 
ttims — m y, to a relative accuracy of < 10 -5 (10~ 6 ) for the 100.pq5+qw3 (200.pq5+qw3) 
slice. Moreover, a quantitative convergence test (cf. appendix E and section XI C; details 
omitted for brevity) shows that these deviations are once again dominated by (the expected) 
4th order finite differencing truncation errors. 

Instead of perturbing the matter variables, we can perturb the initial slice's geometry 
(3- metric and/or extrinsic curvature components). 15 Since our initial data Ansatze always 
take the initial analytical black hole slice to be vacuum, and the York projection oper- 
ator (4) and (21) obviously preserves this, the resulting slice will therefore also always 
be vacuum. We have previously shown (Thornburg (1993), appendix 5) that in the non- 
spherically-symmetric case, suitable geometry perturbations yield (vacuum) dynamic initial 
slices containing black holes surrounded by gravitational radiation. 

However, in spherical symmetry any vacuum slice is necessarily some slice in (a) 
Schwarzschild spacetime. Although a "warped Schwarzschild" slice of this type has no 
true dynamics, it can still be useful as a numerical-relativity test-bed problem: None of the 
test problems suggested by Centrella et al. (1986) are applicable to the spherically sym- 
metric scalar field system, but York (1989) has suggested testing the numerical evolution 
of "hand-warped" slices in Minkowski spacetime. In the spirit of this latter suggestion, our 
"warped Schwarzschild" slices are of interest as alternative test cases. (As will be seen, they 
also offer a particularly good environment for testing an initial data problem algorithm and 
computer code.) 

To assess the accuracy of these slices, we use wims and R a b c dR abcd as diagnostics: mus 
should be constant across each slice, and R a b c dR abcd should be everywhere the same as in 
a matching-mass Schwarzschild spacetime (i.e. the mass-relative I = R a b c dR abcd should be 
unity across the slice). 

For example, consider the daw5c slice, which uses a perturbation applied to A = g rr , 

A ^ A + O.l x Gaussian(r init =20, tr=5) (30) 

This amplitude gives roughly the same warping of the slice (e.g. deviation of Kj i^schw from 
unity) as for the pqw5 slice. 

Figure 15 shows the relative deviations from unity (errors) of mMs/^totai and 
V-^Schw(mtotai) (i- e - the relative errors in ttims and I) across the daw5c slice, computed with 
two different grid resolutions, where m tota i is the total mass of the slice (since the slice is 
vacuum, this is equal to the Misner-Sharp mass within the inner grid boundary). Both 
diagnostics show only very small deviations across the slices, with m = m total to an relative 
accuracy of < 1CT 5 (5xlCT 7 ) for the 100.daw5 (200.daw5) slice, and / = /schw(m total ) to a 
relative accuracy of ^3xl0 -4 (5xl0~ 5 ) respectively. Moreover, a quantitative convergence 
test (cf. appendix E and section XI C; details again omitted for brevity) shows that these 
errors are once again dominated by (the expected) 4th order finite differencing truncation 
errors. 

Finally, note that the computation of I = R abcd R abcd via (14) and (B23)-(B29) depends 
nontrivially on all the state variables (3-metric, extrinsic curvature, and scalar field compo- 



Or we could perturb both matter and geometry. 
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nents) in our computational scheme. The observation that / = ischw(m total ) to high accuracy 
for the vacuum slices, and more generally that / = ischw(m(r)) f°r the vacuum regions of all 
our slices, thus gives a strong test of the overall correctness of a large part (essentially the 
geometry terms in all the equations) of our code, including freedom from most errors in 
either deriving our continuum equations, or programming the numerical computations. 

XII. CONCLUSIONS AND DIRECTIONS FOR FURTHER RESEARCH 

Like many other researchers, we find York's conformal-decomposition technique to offer 
an elegant, reliable, and fairly efficient means of constructing 3 + 1 initial data. In this 
paper we consider the general case of the York decomposition, where (as is required by most 
black-hole-exclusion slicings) K is nonzero and spatially variable throughout most or all of 
the slices, and hence the full nonlinear 4-vector York equations must be solved numerically. 

The main focus of this paper is on the combination of this numerical solution, with a 
set of Ansatze to provide the York decomposition's inputs and inner boundary conditions: 
We begin with a known black hole initial data slice (Ansatz), apply an arbitrary (generally 
constraint- violating) perturbation (Ansatz) to it, then use the York decomposition to project 
the perturbed field variables back into the constraint hypersurface, using a further Ansatz to 
supply inner boundary conditions for the York equations. (Our algorithm also incorporates 
a final numerical coordinate transformation to (e.g.) restore an areal radial coordinate, but 
this is only for convenience.) 

In comparison to other initial data algorithms, the key advantage of this algorithm is its 
flexibility: It's easily applicable to almost any slicing, and so permits the slicing to be chosen 
for convenience in (say) a black-hole-excluded time evolution, rather than being restricted 
by the initial data solver. This algorithm is also equally applicable in spherical symmetry, 
axisymmetry, or in fully 3-dimensional spacetimes, and in both vacuum and non-vacuum 
spacetimes. 

Collectively, this algorithm's Ansatze provide a great deal of flexibility in controlling 
the physical content of the resulting initial data. An interesting area for further research 
would be in trying to better understand the relationship between this physical content, and 
the choices made for the algorithm's various Ansatze. That is, at present we can't a priori 
(i.e. before solving the York equations) predict very much about the physical content of the 
resulting initial data obtained from a given set of Ansatze choices, so if some particular type 
of initial data is desired, some trial and error may be needed to obtain the desired results. A 
better understanding of this would both be of practical benefit in more easily constructing 
initial data with desired properties, and quite likely prove informative in mapping the range 
and phenomenology of possible initial data obtainable using our methods. 

Turning now to the individual Ansatze, the choice of the initial known black hole slice 
seems fairly straightforward: Given the black-hole-exclusion requirement that the slicing 
and all the 3 + 1 field variables be nonsingular near to and on the horizon, the natural choice 
is an Eddington-Finkelstein slice of Schwarzschild spacetime in spherical symmetry, or a 
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Kerr slice of Kerr spacetime in more general cases. 16 

In contrast, the choices for our other Ansatze (the perturbation to apply to the initial 
slice, and the inner boundary conditions for the York decomposition) are less clear-cut. The 
choices presented here (respectively adding a Gaussian to one or a few of the 3-metric, extrin- 
sic curvature, or matter field variable components; and Dirichlet inner boundary conditions) 
are essentially the first ones we tried, and as yet we have made only limited experiments 
with other choices. The choices presented here seem to work adequately, in the sense of 
yielding dynamic black hole spacetimes suitable for numerical evolutions (Thornburg (1993, 
1998)), but clearly these choices are only a tiny subset of the possible parameter space. Both 
numerical and theoretical studies of a broader range of choices would be useful here. 

As well as a general discussion of the initial data problem algorithm on generic slices, in 
this paper we present numerical results for a specific model problem system, the spherically 
symmetric scalar field. We find that the resulting full nonlinear 4-vector York equations can 
be solved robustly and with high accuracy using a Newton-Kantorovich "outer linearization" , 
followed by a standard direct solution of the resulting sequence of "inner" linear elliptic 
PDEs, using 4th order finite differencing on a smoothly nonuniform grid. It seems likely 
that a wide range of other numerical methods would also work well on this problem. 

We haven't encountered any difficulties with the York equations failing to have valid and 
physically reasonable solutions, even with strong perturbations yielding very massive scalar 
field shells (up to ~ 17 times the black hole's mass in our tests). However, we also don't have 
any rigorous proofs that such difficulties couldn't arise for different choices of the various 
Ansatze. Further research in this direction would be valuable. 

As examples of our initial data solver, we discuss a number of Eddington-Finkelstein- 
like initial data slices containing black holes surrounded by scalar field shells, and also 
several vacuum black hole slices with nontrivial slicings. As described in table I, all of our 
test slices use (local) Gaussian perturbations added to one or two of the scalar field, 3- 
metric, or extrinsic curvature components. The slices vary in their perturbations' positions, 
amplitudes, and widths, in which field variable or variables is/are perturbed, and in the 
numerical parameters (grid resolution and outer boundary position). 

We use a number of diagnostics to analyze these slices, including apparent horizon po- 
sitions, the Misner-Sharp and integrated-scalar-field mass functions m MS and m M (respec- 
tively), the scalar field's 3-energy density p and radial 3-energy density 47igggp, the extrinsic 
curvature scalar K, the 3-Ricci scalar R and quadratic curvature invariant RijR % \ the 4-Ricci 
scalar ^R, and the 4-Riemann quadratic curvature invariant R a b cc iR abcd ■ 

The focus of our analysis of the test slices is on exploring the general phenomenology of 
the initial data slices, and on testing and validating the accuracy of the initial data algorithm 
itself. The computed slices' scalar field shells have masses ranging from as low as 0.17 to as 
high as 17 times the black hole mass m^, i.e. the scalar field shells comprise from as little as 
14% to as much as 94% of the total slice masses. The scalar field shells have approximately 



16 Another option, which we haven't explored, would be to use an initial slice constructed by some 
numerical method, either this one or another. Or one could use a non-black-hole initial slice (e.g. a 
slice in Minkowski spacetime) , and try to choose the remaining Ansatze so as to produce a "geon" 
(cf. Abrahams and Evans (1992)). 



24 



Gaussian profiles with fractional widths o"/r center (for the radial scalar field density 4tt geep) 
ranging from 3.3% to 14%, and positions (radii) ranging from 13mt,h to 28mbh- 

In all cases we find that the computed slices are very accurate: For grids with resolutions 
of Ar/r 0.02 (0.01) near the perturbations, the numerically computed energy constraints 
are ^ 10~ 8 (10~ 9 ) in magnitude; the numerically computed momentum constraints are typ- 
ically ~ 3 orders of magnitude smaller. The Misner-Sharp and integrated-scalar-field mass 
functions agree to within relative accuracies of ^ 10~ 6 (10 -7 ) respectively. In tests of nu- 
merically constructed vacuum slices with nontrivial slicings, where we would expect the 
Misner-Sharp mass function to be constant and the 4-Riemann quadratic curvature invari- 
ant R abcd R abcd to be equal to its Schwarzschild-spacetime value, we find these properties 
to hold to within relative errors of ^$ 10~ 5 (5xl0~ 7 ) and <3xl0~ 4 (5xl0~ 5 ) respectively. 
Except for floating point roundoff noise (most prominently in the 3-Ricci scalar R), the 
errors we do see in these and other diagnostics are all generally dominated by the expected 
0((Ar) 4 ) finite differencing truncation errors. 

As well as our main focus on the initial data problem, in the appendices of this paper we 
also present results in several related peripheral areas, including a tabulation of some of the 
key 3 + 1 field variables for an Eddington-Finkelstein slice of Schwarzschild spacetime, the 
formulation of the Misner-Sharp mass function without restrictions on the slicing or spatial 
coordinates, numerical 3-coordinate transformations, and tests of the convergence of finite 
differencing computations to the continuum limit. We also review some little-known, but 
surprising, numerical analysis results on the non-smoothness of the errors incurred when in- 
terpolating data from one grid to another by the usual moving-local-interpolation schemes, 
and discuss these results' import for numerical coordinate transformations and horizon find- 
ing. 
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APPENDIX A: SCHWARZSCHILD SPACETIME IN 
EDDINGTON-FINKELSTEIN COORDINATES 

Although the Eddington-Finkelstein coordinates (t, r, 6, <fi) for Schwarzschild spacetime 
are well known (Misner et al. (1973), box 31.2), so far as we know the 3 + 1 field variables for 
an Eddington-Finkelstein slice of Schwarzschild spacetime haven't been published before in 
the open literature. For reference we list the key ones here; we have given a more complete 
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listing elsewhere (Thornburg (1993), appendix 2). For a Schwarzschild spacetime of mass 



m, 



so that 



a = 7rr? (A1) 
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i^i^ = — ^ Z - 3/ 2 (All) 



Independent of the slicing, Schwarzschild spacetime has p = ^ 4 ^i? = (vacuum) and 
RabcdR abcd = 48m 2 /r 6 (Misner et al. (1973), section 31.7). 

APPENDIX B: EQUATIONS FOR THE SPHERICALLY SYMMETRIC SCALAR 

FIELD SYSTEM 

In this appendix we tabulate the equations for all the 3 + 1 variables involved in our 
initial data algorithm, in terms of the state variables A, B, X, Y, P, and Q. 

1. Geometric Variables 

The only nonzero (3-)Christoffel symbols are 
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The 3-Ricci tensor, scalar, and quadratic curvature invariant are 



Rij 


diag R rr Roe Rgg sin 2 9 






d rr B .(drB) 2 ^drA^drB) 
B 2 B 2 2 




Ree = 


2 A 4 A 2 




R = 


o <9 rrJ B ^drBf (d r A)(d r B) 
AB 2 AB 2 A 2 B 


2 
B 


RijR i:) = 


3 (d rr B) 2 n d rr B 3 (d r A)(d r B)(d„ 
2 A 2 B 2 " AB 2 2 A 3 B 2 


B) 



(Bla) 

(Bib) 

(Blc) 

(Bid) 
(Ble) 
(Blf) 

(Big) 

(B2) 
(B3) 
(B4) 
(B5) 



(d r B) 2 (d„.B) 
A 2 B 3 



2 (d r A)(d r B) 3 (d r A) 2 (d r B) 2 l idrAKdrBf ^Bf 
B 2 A 2 B 2 8 A 4 B 2 2 A 3 B 3 4 A 2 B^ 1 ' 

In spherical symmetry it's relatively easy to locate apparent horizons, by simply zero- 
finding on the apparent horizon equation (York (1989)) 

H = -^S- - 2- = (B7) 

[In detail, to locate apparent horizons we first search the grid function H for sign changes. 
For each sign change (distinct apparent horizon), we construct a local (Lagrange) interpo- 
lating polynomial for the H values in the manner of appendix F, then use the ZEROIN 
code of Forsythe et al. (1977), chapter 7, to locate the interpolating polynomial's zero.] 

Note that in spherical symmetry any apparent horizon must have r area i = 2m. [To 
see this, 17 observe that in spherical symmetry, the Hawking mass function agrees with the 
Misner-Sharp one, and is given by m H = ^r areal (l — ^uj + ujJ), where u + and uj_ are the null 
expansions, normalized to +2 in flat spacetime. An apparent horizon is defined by one or 
both of the expansions vanishing, so m = m H = |^ ar eai there.] 



17 This argument was suggested to the author by N. 6 Murchadha. 
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2. Scalar Field Variables 



The 3 + 1 energy and momentum densities are 



4irf 



-PQ 
PQ 
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and the spatial stress-energy tensor and its trace are 

AnTij = diag [ 4vrT rr ^T ee 4kTm sin 2 9 



AnT„ = \P 2 + \AQ 2 
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(B8) 
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3. Constraints, Mass Functions, and 4-Invariants 



The energy and momentum constraints are 
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In spherical symmetry the mass contained within a given radius is given by the Misner- 
Sharp mass function (Misner and Sharp (1964); Misner et al. (1973), section 23.5). As 
discussed in appendix C, for our purposes this is most usefully expressed in the form 



, r-( Ad r Bf Y 2 ' 



(B18) 



Notice that since B = ggg = 0(r 2 ), the leading term in (B18) is 0(r) for large r. For 
our slices wms is always 0(1) outside the scalar field shells, so the remaining terms in (B18) 
(must) nearly cancel the leading term. This cancellation, and thus our computation of Toms, 
is thus quite sensitive to small finite differencing and roundoff errors at large r. 

As discussed in section VIII, we thus use the integrated-scalar-field mass function m M 
defined by (13). Guven and O Murchadha (1995) define \i as 
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V = \ —Qf- ( r areal^ J J 

where £ is the proper radial distance in the slice, i.e. d£ = J~g^-dr. We thus have 
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To actually calculate m M , we reformulate the integral (13) in terms of our nonuniform 
gridding coordinate w (cf. section VIII), 
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then use the alternate form of Simpson's rule given by Press et al. (1986), (4.1.14), to do the 
numerical integration. (The standard form of Simpson's rule is inconvenient here because it 
only yields results at every 2nd grid point.) 

As noted in section VIII, the 4-Ricci scalar is given by = 8ir(p — T), so we have 
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We compute the 4-Riemann curvature invariant R a bcdR abcd via the general expression (14), 
using the subexpressions 
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The only nonzero components of ViKj k and V\iK^ k are 
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4. York Decomposition Variables and Boundary Conditions 

The terms and coefficients in the York equation (5) and its Jacobian (8) are 
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In spherical symmetry, the outer boundary conditions (10a) and (11) become 
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d r Vl + - = O ( \ ) « (at r = r max ) (B44) 



respectively. 



APPENDIX C: THE MISNER- SHARP MASS FUNCTION IN SPHERICAL 

SYMMETRY 

One of our major diagnostics is the Misner-Sharp mass function, which we use in the 
form (B18). Here we outline the derivation of this result. 

Misner and Sharp (1964) gave the original derivation of the mass function for a spherically 
symmetric dynamic spacetime. However, their results [their equations (6.7) and (6.9)] aren't 
directly usable in a 3 + 1 code such as ours, for two reasons: First, their results are expressed 
in terms of Lagrangian-coordinate (comoving) matter variables, whereas we use Eulerian 
(non-comoving) coordinates. Second, their results are in the form of integrals from the 
origin out to a given radius (analogously to (13) but with lower limits of r = 0), whereas 
with the black hole exclusion technique the field variables aren't known (and in fact are 
generally singular) in the neighbourhood of the origin. 

We thus begin instead with the (equivalent) alternate formulation of the Misner-Sharp 
mass function given by Misner et al. (1973), section 23.5, 

rrius = \f ( 1 - — ) (CI) 

^ \ 9rr / 

where the bars denote the use of Schwarzschild coordinates x a = (i,f,0,(f>), defined by 
requiring the 4-metric to be diagonal (gj f = 0) and the radial coordinate to be areal (gee = 
f 2 ). 

In our 3 + 1 code we permit the coordinates to be arbitrarily chosen (within the re- 
striction of spherical symmetry), i.e. we permit the 4-metric to be non-diagonal and/or the 
radial coordinate to be non-areal. To transform (CI) back to such generic coordinates, we 
first define an intermediate coordinate system x a = (t,f,9,(f>), which uses our generic time 
coordinate t but retains the areal radial coordinate f. From the requirement that gif — 0, a 
straightforward calculation shows that g ff = g-w — (9tf) 2 /gtt, and hence 

1_ 

m M s = 2 r 

[This result appears somewhat problematic, as it has m apparently depending on the freely 
specifyable ^-coordinate lapse and shift a and f3. However, the areal radial coordinate 
condition gee = r 2 actually fixes the ratio [3 /a which enters into m (i.e. this ratio actually 
isn't freely specifyable), so there's no difficulty] 

From the requirement that gee = r 2 it then follows that dr/df = 2\fB /d r B, and hence 

-ms = \Vb(i- + (C3) 

which is the desired generic-coordinates form for the mass function. 
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APPENDIX D: NUMERICAL COORDINATE TRANSFORMATIONS 



In 3 + 1 numerical relativity one often uses coordinate conditions which fix the form of 
the spatial metric. For example, one may wish to force the radial coordinate to be areal 
{dee — r ' 2 )- Even if our Ansatz for the base field variables satisfies such a condition, in 
general the York projection operator (4) doesn't preserve it. As the final step in our overall 
initial data algorithm (step 4 in our presentation of section VII), we thus optionally apply 
an explicit numerical coordinate transformation to the output of the York projection, so as 
to (e.g.) convert it back to an areal radial coordinate. 

We actually use a numerical-coordinate-transformation algorithm somewhat more gen- 
eral than is needed just for this purpose: our algorithm actually applies a generic (nonsin- 
gular and invertible) numerically specified radial-coordinate transformation r — > f = f(r). 
(To specify a transformation to an areal radial coordinate, we simply take f = y/gee.) 

Numerical coordinate transformations are conceptually straightforward, but somewhat 
complicated in detail due to the field variables only being known on the finite difference 
grid. In our case the use of the nonuniform gridding coordinate w further complicates 
matters. Figure 16 summarizes our numerical coordinate transformation algorithm; given 
this flowchart the operations involved should be self-evident. Notice that as well as tensor- 
transforming the field variables, the algorithm must also interpolate them from the r to the 
f grid. We use moving local (Lagrange) polynomial interpolation for this, as discussed in 
detail in appendix F. 

APPENDIX E: CONVERGENCE TESTS 

As discussed in section XI C, several different error sources are present in a 3 + 1 code 
such as ours, but the dominant one is usually finite differencing truncation error. As has 
been forcefully emphasized by Choptuik (1986, 1991); Choptuik et al. (1992), in such a 
situation a careful comparison of the errors in approximating the same physical system at 
different grid resolutions can yield a great deal of information about, and very stringent 
tests of, the code's numerical performance and correctness. 

In the present context, consider some diagnostic Z whose true (continuum) value Z* is 
known. [For example, the true values of the energy and momentum constraints C and C l 
are zero, and in section XI E the true value of ///schw(m total ) in the (vacuum) daw5 slice is 1, 
where I = R a t, c dR abcd ■] Suppose we have a pair of numerically computed approximate initial 
data slices, identical except for having a 2:1 ratio of grid spacings Aw. As discussed in detail 
by Choptuik (1991), if all the field variables are smooth and the code's numerical errors are 
dominated by truncation errors from nth order finite differencing, the numerically computed 
diagnostics Z must satisfy the Richardson expansion 

Z[Aw] =Z* + {Aw) n f + 0((Aw) n+2 ) (Ela) 
Z[Aw/2] = Z* + (Aw/2) n f + 0((Aw) n+2 ) (Elb) 

at each grid point, where Z[Aw] denotes the numerically computed Z using the grid spacing 
Aw, and / is an 0(1) smooth function depending on various high order derivatives of Z* 
and the field variables, but not on the grid resolution. [We're assuming centered finite 
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differencing here in writing the higher order terms as 0((Aw) n+2 ), otherwise they would 
only be 0((Aw) n+1 ).] Neglecting the higher order terms, i.e. in the limit of small Aw, 
we can eliminate / to obtain a direct relationship between the code's errors at the two 
resolutions, 

Z[Aw/2] - Z* = 1_ 
Z[Aw] - Z* 2 n { ' 

which must hold at each grid point common to the two grids. 

We use two methods for testing how well or poorly any particular set of (finite-resolution) 
numerical results satisfies the convergence criterion (E2): One method, used here in figures 8, 
11, 14, and 15, is to plot the absolute values of the low-resolution and high-resolution errors 
(|Z[Aw] — Z*\ and |Z[Aiu/2] — Z*\ respectively) on a common logarithmic scale. If, and 
given the arguments of Choptuik (1991), in practice only if, the error expansions (El) are 
valid with the higher order error terms negligible, i.e. if and only if the errors are indeed 
dominated by the expected nth order finite difference truncation errors, then the two curves 
will be identical except for the high-resolution curve being offset downwards (towards smaller 
errors) by a factor of 2™ from the low-resolution curve. 

Our other method (Thornburg (1993, 1996)) for testing our results against the conver- 
gence criterion (E2), used here in figure 9, is to plot a scatterplot of the high-resolution errors 
Z[Aw/2] — Z* against the low-resolution errors Z[Ato] — Z* at the grid points common to 
the two grids. If, and again in practice only if, the error expansions (El) are valid with the 
higher order error terms negligible, i.e. if and only if the errors are indeed dominated by the 
expected nth order finite difference truncation errors, then all the points in the scatterplot 
will fall on a line through the origin with slope 1/2™ (= ^ for 4th order convergence). 

Comparing these analyses, each has its advantages and disadvantages: The logarithmic 
analysis simultaneously tests the convergence criterion (E2), and gives a good qualitative 
overview of the overall magnitude and spatial variation of the errors. Moreover, the use of a 
logarithmic plot gives this method excellent dynamic range, i.e. it can easily test deviations 
from the convergence criterion in parts of the grid where the errors are small, even if elsewhere 
in the grid the errors are large. However, this analysis is relatively insensitive to small 
deviations from the expected convergence factor. In contrast, the scatterplot analysis has 
excellent sensitivity to small deviations from the expected convergence factor, but it doesn't 
give any direct indication of where in the grid any problems lie, and it has limited dynamic 
range unless multiple plots at different scales are used (as in figure 9). 

It's important to note that both analyses are applied pointwise, i.e. independently at each 
grid point common to the two grids. This makes these analyses much more sensitive than 
a simple comparison of gridwise norms, which would tend to "wash out" any convergence 
problems occurring only in a small subset of the grid points (e.g. near a boundary). 

[One "false alarm" problem of this type to which the scatterplot analysis is particularly 
subject, is that of boundary finite differencing effects: Near a grid boundary, we use off- 
centered molecules for finite differencing which, while still 4th order in accuracy, are an 
0(1) factor less accurate than centered molecules we use in the grid interior. This in itself 
doesn't cause a problem, but the choice of centered vs. off-centered molecules is based on the 
distance in grid points from the grid boundary, and so may differ at a near-boundary grid 
point between different-resolution grids. For example, with a 2:1 ratio of grid resolutions, a 
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point which is 2 grid points from the boundary (and thus can use a centered 5 point molecule) 
in the high-resolution grid, is only 1 grid point from the boundary (and thus must use an 
off-centered molecule) in the low-resolution grid. The different 0(1) factors in the different 
molecules' truncation errors then cause spurious apparent changes in the convergence ratios 
in (E2) at such grid points. This is the cause of the boundary artifacts in figures 8, 11, 14, 
and 15, and of the outlier points in figure 9.] 

Finally, note that the parameter n, the order of the convergence, is (or at least should 
be) known in advance from the form of the finite differencing scheme. Thus the factor 
of 2 n convergence ratio in (E2), manifested in the vertical offset between the two plots 
in the logarithmic analysis, or in the slope of the line through the origin in the scatterplot 
analysis, isn't fitted to the data points, but is rather an a priori prediction with no adjustable 
parameters. A convergence test of either type is thus a very strong test of the validity of 
the overall finite differencing scheme and the error expansion (El). 

APPENDIX F: THE NON-SMOOTHNESS OF INTERPOLATION ERRORS 

As discussed in appendix D, interpolation plays an important role in numerical coordi- 
nate transformation, and it's also important in horizon finding (Baumgarte et al. (1996); 
Thornburg (1996); Anninos et al. (1996); Huq (1996)). In this appendix we discuss some 
surprising properties of the errors in interpolation. For pedagogical convenience we focus 
only on (Lagrange) polynomial interpolation, but our results actually hold for a wide range 
of interpolation schemes. 

To understand the behavior of interpolation errors, it's useful to consider an artificial 
example where the data being interpolated is in fact already known analytically: Suppose 
we have a smooth continuum function / and a set of grid points {xk} in its domain. To 
simplify our discussion we take the grid points {xk} to be uniformly spaced with spacing 
Ax, but again this restriction isn't essential. 

Now define an interpolating function I using the function values {f(x k )} in the usual 
manner via an n + 1 point moving local interpolant: On each interval [xj,Xj + i] between 
adjacent grid points, we set / to the (unique) nth-degree interpolating polynomial P[i] 
which agrees with the {/(x*;)} values at the n + 1 nearest-neighbour grid points located 
symmetrically about the interval [xj, Xj+i]. 

Consider the interpolation error e = / — / in taking I as an approximation to /, and 
this error's behavior as the grid is refined (Ax — > 0). It's well known that for sufficiently 
smooth /, e = 0((Ax) n+1 ), so the error can be made small by refining the grid. But the 
error function e = e(x) is not smooth! In fact, even if / is analytic, the error function has a 
discontinuous (1st) derivative at each grid point. This holds regardless of the interpolation 
order n. 

The reason for this behavior is that although in any grid interval [xk,Xk+i], I = P[k] 
is a smooth function (an nth-degree polynomial), / is a different such polynomial in each 
separate interval [xk,Xk+i]. At each grid point the two adjacent interval's interpolating 
polynomials agree in their values (so I is continuous there), but not in their derivatives 
(so dl/dx and hence de/dx has a jump discontinuity there). In other words, e is a "bump 
function": it vanishes at each grid point, but between each pair of grid points it rises to a 
nonzero value, here with amplitude 0((Ax) n+1 ). 
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Now suppose we approximate derivatives of / by derivatives of the interpolating function 
/. Because of e's non-smoothness, the Richardson-type argument outlined by Choptuik 
(1991) doesn't apply here. Instead, by considering the behavior of the interpolation process 
for polynomials / of varying degrees, it's straightforward to show that for the pth derivative 
the error in this approximation is 

dP 1 ^dH_<Pl = 0{{Ax) n+l- P) (pi) 

dxP dxP dxP u ' ' v ' 

Hence provided p < n, the error in approximating derivatives this way can still be made 
small by refining the grid, but due to the non-smoothness of the error function the effective 
order of convergence is lowered by the number of derivatives taken. 

This analysis assumes that all differentiation (e.g. to compute d p I/dx p ) is done exactly 
(analytically). However, if local numerical differentiation (finite differencing) is used instead, 
then it's easy to show that the error expansion (Fl) still holds. 

As an example of interpolation errors, consider the function / defined by f(x) = 
exp(sin(|x)), and take n = 3, so the local interpolation uses cubic polynomials fitted to 
4 nearest-neighbour function values. Figure 17(a) shows the error e in this interpolation 
for grids with spacings in the 2:1 ratio Ax = 0.1:0.05, plotted for the range [0,0.5]. (In 
each case the grids are enough larger than this range so the interpolations never encounter 
the grid endpoints.) The overall bump-function shape of the error is clearly evident, with 
e vanishing at the grid points (the interpolation is exact there), and rising to a maximum 
between each pair of grid points. The amplitudes of the "bumps" for the two grid spac- 
ings are in approximately a 16:1 ratio, confirming that we have 4th order convergence here, 
e = 0((A:r) 4 ). 

Figure 17(b) shows the error de/dx in taking dl/dx as an approximation to df/dx for 
this same example, de/dx has a jump discontinuity at each grid point, corresponding to 
the abrupt change in slope going from one "bump" to the next in e (figure 17(a)). The 
amplitudes of de/dx for the two grid spacings are in approximately an 8:1 ratio, confirming 
that we only have 3rd order convergence here, e = 0((Ax) 3 ). 

Our conclusion from this argument is that by virtue of (Fl), if interpolated data is to be 
used as input to further finite differencing computations, then either the interpolation order 
should be raised sufficiently to compensate for the non-smoothness, or Hermite, spline, or 
other explicitly smoothness-preserving interpolation schemes should be used. 

In our present computational scheme we use moving local polynomial interpolation in 
the manner described above. Our overall computational scheme is intended to be 4th order 
accurate, and we use numerical differentiation up to 2nd order (e.g. in the 3-Ricci tensor), so 
we use 5th order (Lagrange) polynomials for local interpolation in our numerical coordinate 
transformation and horizon finding algorithms. These interpolants give 6th order local 
truncation errors, so 2nd numerical derivatives of the interpolation results should still be 
4th order accurate. 
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FIG. 1. This figure shows the behavior of our nonuniform gridding coordinate w. Part (a) shows 
the grid spacing Ar for grids with resolutions of (top to bottom) Aw = 0.01, 0.005, and 0.0025. 
The diagonal dashed lines labeled along the top and right of the figure show lines of constant 
relative grid spacing Ar/r. Part (b) shows the actual w coordinate. In both parts of the figure the 
vertical dashed lines show the outer grid boundaries for w mSLX = 4, 10, and 30. 

FIG. 2. This figure shows the York-transformation conformal factor and vector potential 
for the 200.pqw5 slice. The vertical dashed line at r- m i t = 2 shows the input horizon position. 

FIG. 3. This fi gure shows the relative 3- metric and extrinsic curvature components A = g r n 
B = ggg, X = K rr , and Y = Kgg, for the 200.pqw5 slice. (Note that 5/Sschw = 1 by the definition 
of an areal radial coordinate.) The vertical dashed line just inside r = 2 shows the horizon position. 

FIG. 4. This figure shows the relative 3-invariants K, R, and J = RijR^, for the 200.pqw5 
slice. The horizontal dashed line shows the "zero axis" K/Kschw = R/Rschw = J/Jschw = 1- The 
vertical dashed line just inside r = 2 shows the horizon position. 

FIG. 5. This figure shows the scalar field and mass distributions for the 200.pqw5 slice. Part (a) 
shows the full data range, with logarithmic scales for most of the variables, while part (b) shows 
only the inner part of the grid, with linear scales for all variables. (In part (b) the left axis uses 
a different scale for |( 4 ).R| than for inBp and \P\). In both parts of the figure the vertical dashed 
line just inside r = 2 shows the horizon position. 

FIG. 6. This figure shows the mass-relative 4-Riemann curvature invariant I = R a b c dR abcd for 
the 200.pqw5 slice. The vertical dashed line just inside r = 2 shows the horizon position. Notice 
that to within the accuracy of the plot, I / 'ischw(m(r)) = 1 everywhere outside the scalar field shell. 

FIG. 7. This figure shows the mass-relative 3-invariants K, R, and J = RijR % \ for the 200.pqw5 
slice. The upper and lower horizontal dashed lines show the "zero axes" K/K Schw ^ m ^ = 1 (left 
scale) and R/ Rs c hw(m(r)) = J I ^Schw (m(r)) = 1 (right scale) respectively. The vertical dashed line 
just inside r = 2 shows the horizon position. 

FIG. 8. This figure shows the magnitudes of the numerically computed energy constraint C 
(part (a)), and momentum constraint C r (part (b)), at various points in our initial data algorithm 
for the 100.pqw5 and 200.pqw5 slices. |Cs c hw| (|Cschwl) ^ s magnitude of the energy (momen- 
tum) constraint of the initial Schwarzschild slice before applying the perturbation (step 1 in our 
overall initial data algorithm of section VII) ; | C pertur b | ( | Cp er t U rb I ) ^ s * ne magnitude of the energy 
(momentum) constraint after applying the perturbation (step 2); and |Cfi na i| (|Cg nal |) is the mag- 
nitude of the final energy (momentum) constraint after applying the York algorithm (step 3) and 
the numerical coordinate transformation back to an areal radial coordinate (step 4) . The scale bars 
show a factor of 16 in \C\ (|C r |), for comparison with the vertical spacing between the 100.pqw5 
and 200.pqw5 curves. For the 200.pqw5 |C pcrtur b| curve, only every 5th grid point is plotted. 

FIG. 9. This figure shows a scatterplot of Cfi na i between the 100.pqw5 and 200.pqw5 slices. 
Part (a) shows all the Cfi na i values, while part (b) shows an expanded view of the central region (the 
grid- interior points) of part (a). In each part the line through the origin has slope ^; As discussed 
in appendix E, for 4th order convergence we expect all the points (except for grid-boundary outliers) 
to fall on the line. 
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FIG. 10. This figure shows various embedding and curvature diagnostics for the outermost part 
of the grid in the 100o30.pqw5 slice (part (a)), and the 200o30.pqw5 slice (part (b)). Each part of 
the figure shows the relative K; the mass-relative J = RijR^ and / = R a b c dR abcd ", the relative R; 
and the relative Rst (an empirical least-squares fit to R) perturbed by simulated roundoff errors 
in the 3-metric components, with e = 6xlCT 19 . One outlier point in part (b), that for R/Rs c hw at 
the outermost grid point, falls outside the plot's range and is omitted. (This point is anomalous 
due to boundary finite differencing effects.) Note that the plot key (the choice of line and point 
types used for the different diagnostics) differs between this and the other figures. 

FIG. 11. This figure shows the magnitude of the numerically computed energy constraint, |C|, 
at various points in our initial data algorithm for the 100.pqw5i and 200.pqw5i slices, analogously to 
figure 8(a) for the 100.pqw5 and 200.pqw5 slices. |Cs c hw| is the magnitude of the energy constraint 
of the initial Schwarzschild slice before applying the perturbation (step 1 in our overall initial data 
algorithm of section VII); | C pcrtur b I is the magnitude of the energy constraint after applying the 
perturbation (step 2), and |Cfi na i| is the magnitude of the final energy constraint after applying 
the York algorithm (step 3) and the numerical coordinate transformation back to an areal radial 
coordinate (step 4). The scale bar shows a factor of 16 in |C|, for comparison with the vertical 
spacing between the 100.pqw5i and 200.pqw5i curves. For the 200.pqw5i |C pe rturb| curve, only 
every 5th grid point is plotted. 

FIG. 12. This figure shows the scalar field and mass distributions for the 200.pqw5b (part (a)) 
and 200.pqw5c (part (b)) slices, analogously to figure 5(a) for the 200.pqw5 slice. In both parts of 
the figure the left scale is identical to that of figure 5(a); the right scale differs for each part. In 
part (a) the vertical dashed line just inside r = 2 shows the horizon position. (As discussed in the 
text, there is no horizon within the numerical grid in the 200.pqw5c slice plotted in part (b).) 

FIG. 13. This figure shows the scalar field and mass distributions for the 400.pqwl slice, anal- 
ogously to figure 5(a) for the 200.pqw5 slice and figure 12 for the 200.pqw5b and 200.pqw5c slices. 
Part (a) shows the full data range, while part (b) shows only the inner part of the grid. In both 
parts of the figure the vertical dashed line just inside r = 2 shows the horizon position. 

FIG. 14. This figure shows the relative deviation \m,Ms/ m n ~ M °^ the Misner-Sharp mass 
function myis from the integrated-scalar-field mass function m^, for the 100.pw5+qw3 and 
200.pw5+qw3 slices. The scale bar shows a factor of 16 in the deviations, for comparison with the 
vertical spacing between the 100.pw5+qw3 and 200.pw5+qw3 curves. The vertical dashed line just 
inside r = 2 shows the horizon position. 

FIG. 15. This figure shows the relative errors |mMs/ m totai — 1| and |-f/ischw(m tota i) — 1| i n the 
Misner-Sharp mass function tbms and the 4-Riemann curvature invariant I = R a i )C( j i R abcd for the 
100.daw5c and 200.daw5c slices, where m to tai is the total mass of the slice. The scale bar shows 
a factor of 16 in the errors, for comparison with the vertical spacing between the 100.daw5c and 
200.daw5c curves. The vertical dashed line just outside r = 2 shows the horizon position. 

FIG. 16. This figure shows the overall data flow in our numerical coordinate transformation 
algorithm. Oval boxes denote grid functions, rectangular boxes denote operations performed on 
grid functions, solid arrows show data flows presently used in our algorithm, and the dashed 
arrow shows an additional data flow which would be needed to transform contravariant tensors in 
addition to covariant ones, f @ w and f @ w denote a grid function / defined on the w and w 
grids respectively. 
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FIG. 17. This figure shows the interpolation error e = I — f (part (a)) and its derivative 
de/dx = dl/dx — df/dx (part (b)), for the function / defined by f(x) = exp(sin(|x)), cubically 
interpolated to grids with spacings Ax = 0.1 and 0.05. The plot key for part (a) also applies to 
part (b). 
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TABLE I. This table gives various parameters for our test initial data slices. For each slice 
the table shows a descriptive name, the spatial grid resolution and size, the initial perturbation 
(step 2 in our overall initial data algorithm of section VII), an approximate Gaussian fit to the 
slice's radial scalar field density 4wBp, the horizon's areal radius, and a summary of the slice's 
mass distribution (broken down into black hole, scalar field, and total mass). 



Name 



Aw 



Ar/r 
at r=20 



w 



max ' max 



Initial Perturbation 



Approximate Pre 
Scalar Field Shell 



100.pqw5 

200.pqw5 

100ol0.pqw5 

200ol0.pqw5 

100o30.pqw5 

200o30.pqw5 

100.pqw5i 

200.pqw5i 

100.pqw5b 

200.pqw5b 

100.pqw5c 

200.pqw5c 

200.pqwl 

400.pqwl 

100.pw5+qw3 

200.pw5+qw3 

100.daw5c 

200.daw5c 



0.01 

0.005 

0.01 

0.005 

0.01 

0.005 

0.01 

0.005 

0.01 

0.005 

0.01 

0.005 

0.005 

0.0025 

0.01 

0.005 

0.01 

0.005 



0.02 

0.01 

0.02 

0.01 

0.02 

0.01 

0.02 

0.01 

0.02 

0.01 

0.02 

0.01 

0.01 

0.005 

0.02 

0.01 

0.02 

0.01 



248 



10 813 



> P->P + 0.02 xGaussian(r imt =20,cr=5) 0.0735 xGaussian(r= 21. 



30 2 780 

4 248 P^P + 0.02 xGaussian(r imt =10,<7=5) 0.0205 xGaussian(r= 12. 

4 248 P^P + 0.05 xGaussian(r mit =20,o-=5) 0.0365 xGaussian(r= 24. 

4 248 P^P + 0.1 xGaussian(r init =20,o-=5) 0.91 xGaussian(r= 31. 



4 248 
4 248 
4 248 



P^P + 0.1 xGaussian(r imt =20,<7=l) 1.42 xGaussian(r=22.' 



' P -> P + 0.02 x Gaussian(r init =20, cr=5) ' 
_Q^Q + 0.03 x Gaussian(r illit =20, cr=3) . 

A^A + 0.1 xGaussian(r ini t=20, cr=5) 



0.225 xGaussian(r= 22. 



none 



a As discussed in the text in section XI E, there's no apparent horizon within this slice's numerical 
grid; the "black hole" mass shown (within parentheses) is actually the mass within the inner grid 
boundary. 

b The initial perturbation is to both P and Q for this slice. 
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TABLE I. This table gives various parameters for our test initial data slices. For each slice 
the table shows a descriptive name, the spatial grid resolution and size, the initial perturbation 
(step 2 in our overall initial data algorithm of section VII), an approximate Gaussian fit to the 
slice's radial scalar field density 4irBp, the horizon's areal radius, and a summary of the slice's 
mass distribution (broken down into black hole, scalar field, and total mass). 



Ar/r 



Name 


A 


w 


at r = 20 


100.pqw5 





01 





02 


200.pqw5 





005 





01 


100ol0.pqw5 





01 





02 


200ol0.pqw5 





005 





01 


100o30.pqw5 





01 





02 


200o30.pqw5 





005 





01 


100.pqw5i 





01 





02 


200.pqw5i 





005 





01 


100.pqw5b 





01 





02 


200.pqw5b 





005 





01 


100.pqw5c 





01 





02 


200.pqw5c 





005 





01 


200.pqwl 





005 





01 


400.pqwl 





0025 





005 


100.pw5+qw3 


01 





02 


200.pw5+qw3 


005 





01 


100.daw5c 





01 





02 


200.daw5c 





005 





01 



w max ''max 



Initial Perturbation 



Approximate Profile of 
Scalar Field Shell (4ttB P ) 



Mass 



Horizon 
Radius T>H + SF~ 



Total 




4 248 



4 248 



4 248 



4 248 



4 248 



4 248 



P->P + 0.02xGaussian(nnit=20, (7 = 5) 0.0735X Gaussian(r= 21.8 , (7= 3.5 ) 1.952 0.976 + 0.641 = 1.617 



P P + 0.02 X Gaussian(n n it = 
P — ► P + 0.05 X Gaussian(n n it : 
P^P + 0.1 X Gaussian(r;nit : 
P^P + 0.1 X Gaussian(r;nit : 



= 10,(7 = 5) 0.0205xGaussian(r= 12.4 ,a= 3.2 ) 1.976 0.988 + 0.167 = 1.155 
= 20,(7 = 5) 0.0365xGaussian(r= 24.4 ,(7= 3.7 ) 1.733 0.866 + 3.405 = 4.271 
= 20,(7 = 5) 0.91 xGaussian(r= 31.0 ,(7= 4.2 ) 



(0.576) + 9.737 = 10.313 



P P + 0.02X Gaussian(r in it =20, (7 = 5) 
Q^Q + 0.03X Gaussian(r in it =20, (7 = 3) 



A^A + 0.1 xGaussian(n n it=20, (7 = 5) 



20,(7 = 1) 1.42 xGaussian(r = 22. 71, (7 = 0.76) 1.784 0.892 + 2.692 = 3.584 

b 

0.225 xGaussian(r= 22.1 ,(7= 2.5 ) 1.872 0.936 + 1.666 = 2.602 

none 2.043 1.022 - 1.022 
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